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1  INTRODUCTION 

1.1  Purpose  and  Scope 

The  purpose  of  this  report  is  to  document  the  software  implementation  of  a  Karhunen- 
Loeve  estimator  of  local  surface  gravity  fields  from  measurements  of  airborne  gravity  gra¬ 
dient  fields,  and  to  present  and  discuss  the  results  of  associated  numerical  experiments. 
Extensive  derivations  and  theoretical  justifications  are  avoided.  For  these  the  reader  is 
referred  to  Bose  [2]  of  which  a  brief  extract  covering  the  principal  equations  is  presented 
in  Appendix  A. 

Here  we  place  special  emphasis  on  the  validation  of  the  software  using  synthetic  data 
from  a  single  dipole  gravitational  field,  which  we  simulated  ourselves.  The  main  body  of 
experiments  were  performed  on  simulated  measurements  produced  at  the  Naval  Surface 
Weapons  Center  (NSWC)  based  on  multiple  layers  of  random  dipoles. 

1.2  Problem  Statement 

Given  a  measurement  of  the  gravity  gradient  field  collected  on  a  regular  grid  at  an  altitude 
h  above  the  earth’s  surface,  estimate  the  gravity  vector  field  at  altitude  0.  Design  and 
implement  in  software  an  estimator  based  on  the  Karhunen-Loeve  algorithm.  Validate  the 
software  using  controlled  data  and  test  it  on  foreign  data  supplied  by  the  NSWC.  Perform 
an  error  analysis  on  i  he  results. 

1.3  Overview  of  Report 

Chapter  2  briefly  discusses  the  technical  approach  including  an  overview  of  the  algorithm, 
software  design,  various  sets  of  measurements,  and  finally  error  analysis.  The  theory  and 
results  of  validating  the  software  using  known  single  dipole  simulations  are  presented  in 
Chapter  3.  Naval  Surface  Weapons  Center  (NSWC)  data  and  related  pre-processing  are 
discussed  in  Chapter  4.  Results  of  testing  the  Karhunen-Loeve  estimator  for  various  sets  of 
measurements  are  presented  in  Chapter  5.  Observations  made  from  analysing  the  results 
are  discussed  in  Chapter  6.  A  summary  of  the  research  performed,  significant  conclusions 
drawn  and  recommendations  for  future  research  are  outlined  in  Chapter  7. 
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2  TECHNICAL  APPROACH 

2.1  Algorithm  Overview 

The  modeling  approach  taken  here  is  to  exploit  the  marnage  of  physical  geodesy  and 
random  process  theory.  Laplace’s  equation  is  solved  with  the  unknown  mass  distribution 
below  the  surface  of  the  earth  modelled  as  a  two-dimensional  white  noise  layer  representing 
the  vertical  derivative  of  the  disturbance  potential  to  any  pre-specified  order.  This  results 
in  a  series  solution  of  the  disturbance  potential  wherein  the  unknown  coefficients  of  the 
expansion  are  forced  to  be  uncorrelated  by  invoking  the  Karhunen-Loeve  condition.  It 
can  be  shown  that  the  disturbance  potential  covariance  obtained  from  this  model  is  both 
non-stationary  and  non-isotropic. 

The  six  (6)  gravity  gradient  measurements  are  represented  in  terms  of  the  Karhunen- 
Loeve  scries  expansion  of  the  disturbance  potential  resulting  in  six  basis  functions.  These 
basis  functions  are  shown  to  be  orthogonal.  A  linear  meap  square  estimator  utilizing 
all  the  gravity  gradient  measurements  simultaneously  is  obtained  in  continuous  domain 
by  solving  an  integral  equation  involving  the  estimator  gains  which  are  represented  by  the 
same  orthogonal  basis  functions.  The  discrete  implementation  of  the  estimator  is  facilitated 
by  exploiting  the  orthogonality  or  near-orthogonality  of  the  transformation  matrices  so 
that  inverting  these  matrices  becomes  computationally  trivial.  It  turns  out  that  the  near- 
orthogonal  cosine  matrix  must  be  of  even  order  for  its  inverse  to  exist,  a  minor  restriction. 
The  gravity  disturbance  vector  is  obtained  in  a  two-dimensional  grid  which  can  be  denser 
than  the  measurement  grid  and  also  at  any  altitude,  including  the  surface  of  the  earth. 
Thus,  interpolation  and  downward  continuation  are  performed  automatically.  Details  of 
the  mathematical  development  of  the  algorithm  are  discussed  in  3ose  [2].  A  summary 
extraction  of  the  relevant  equations  necessary  for  software  development  are  outlined  in 
Appendix  A. 

2.2  Software  Overview 

A  modular  software  package  was  constructed  to  realize  the  estimation  algorithm  and  facil¬ 
itate  numerical  experiments,  see  Figure  1.  The  estimation  process  is  split  into  two  parts: 
in  module  KLE  (Karhunen-Loeve  Estimator)  the  measured  gravity  gradient  fields  are  pro¬ 
cessed  to  produce  the  K-L  (Karhunen-Loeve)  coefficients  &ki  of  the  estimated  potential 
field;  in  module  SYN  (synthesizer)  the  coefficients  .  re  used  to  synthesize  the  gravity  vec¬ 
tor  field  from  the  now  known  coefficients.  Incidentally,  the  SYN  module  can  be  used  to 
simulate  random  gravitational  fields  by  feeding  it  coefficients  from  a  suitably  distributed 
psuedo-random  variable.  These,  two  modules  are  constructed  so  that,  in  principle  at  least, 
any  combination  of  derivatives  of  the  potential  may  be  used  as  the  observation,  and  any 
combination  may  be  estimated:  this  includes  the  potential  T,  the  gravity  vector  T,  and  the 

gradient  tensor  T  or  any  subset  of  their  components.  For  computational  speed,  fast  sine 
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and  cosine  transforms  (FST  and  FCT)  are  used  to  compute  the  discrete  sine  and  cosine 
transforms  (DST  and  DCT)  in  both  the  KLE  and  SYN  modules.  Note  that  the  sine  matrix 
of  our  algorithm  is  identical  to  the  standard  DST,  whereas  our  cosine  matrix  differs  from 
the  standard  DCT  in  that  it  lacks  the  first  and  last  rows  and  columns  corresponding  to 
the  region  boundaries.  This  fact  does  not  prevent  the  use  of  the  FST  and  FCT  to  obtain 
the  speed  advantage  of  the  fast  transforms.  In  order  to  fully  realize  this  speed  advantage, 
the  number  of  samples  transformed  plus  one  must  be  .  .  (odd)  product  of  small  primes. 

In  addition,  th'ure  are  modules  which  prepare  the  measurements  for  processing  and 
which  generate  plotting  files  from  any  of  the  measured  or  computed  fields.  Yet  another 
module  was  created  to  merge  the  K-L  coefficients  from  west-east  and  north-south  tracks  to 
form  K-L  coefficients  representing  the  best  estimate  from  all  of  the  data,  i.e.  the  S-N  and 
W-E  grids  together.  A  module  was  designed  to  generate  single  dipole  gravitational  fields 
including  the  gravity  vector  and  gradients;  the  controlled  data  produced  by  this  module 
were  used  to  exercise  and  validate  the  estimation  software  during  its  development. 

These  modules  are  coded  in  portable  Fortran  77  and  are  invoked  by  VAX  (tm)  DCL 
command  procedures  to  do  production  rims.  Internally,  all  the  programs  compute  using 
physically  consistent  SI/MKS  units.  Inputs  and  outputs  may  be  expressed  in  units  con¬ 
venient  for  the  user,  e.g.  mgal  or  Eotvos.  Any  conversions  between  internal  and  external 
units  are  performed  when  the  data  are  input  or  output. 


Figure  1:  Software  Modules 
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2.3  Single  Channel  Measurements  on  a  Single  Grid 

Experiments  were  performed  consisting  of  estimating  the  gravity  vector  from  each  gravity 
gradient  component.  The  measurement  data  were  presented  on  a  rectangular  grid.  File  1 
measurements  axe  taken  every  1  km  on  south-north  (S-N)  tracks  which  are  5  km  apart. 
File  2  measurements  were  taken  on  west-east  (W-E)  tracks  with  the  same  spacings. 

Figures  2  and  3  illustrate  rectangular  grids  with  S-N  and  W-E  tracks,  respectively.  The 
sample  points  used  in  the  estimation  are  represented  by  the  drts.  The  border,  which  is 
not  processed,  is  represented  by  the  enclosing  box.  The  width  and  height  of  the  box  are 
the  half  periods  A  and  B  of  the  Fourier  series  representation  of  the  fields.  The  numbers 
of  interior  sample  points  horizontally  and  vertically  are  given  by  K  and  L.  The  numbers 
of  samples  in  the  illustrations  are  not  those  of  the  actual  data  processed. 


Figure  2:  Rectangular  Grid:  File  1,  S-N  Tracks 


Figure  3:  Rectangular  Grid:  File  2,  W-E  Tracks 


4 


2.4  Full  Tensor  Measurement  on  a  Single  Grid 

After  successful  estimation  of  the  gravity  vector  from  single  gravity  gradient  components, 
the  full  gravity  gradient  tensor  was  used  as  a  measurement  and  an  estimate  successfully 
produced.  The  six  components  measured  were  Tir,  Tyy,  Tzz,  Tiy,  Tyz,  and  Txz.  Strictly 
speaking,  any  one  of  the  first  three  of  these  may  be  considered  redundant  in  view  of 
Laplace’s  equation,  but  all  were  still  included  and  all  six  measurement  noises  were  consid¬ 
ered  independent.  For  simplicity,  all  of  the  six  measurement  noise  variances  were  taken  to 
be  equal  to  unity.  This  experiment  was  performed  with  both  file  1  and  file  2  data,  i.e.  the 
same  measurements  which  were  used  singly  in  the  previous  set  of  experiments. 

2.5  Full  Tensor  Measurement  on  a  Dual  Grid 

Files  1  and  2  contain  measurements  of  the  same  fields  over  the  same  survey  region,  but 
one  has  south-north  tracks  and  the  other  has  east-west  tracks.  Samples  are  taken  1  km 
apart  and  the  tracks  are  5  km  apart.  Thus  the  sampling  densities  differ  in  each  horizontal 
direction.  We  would  like  to  use  both  sets  of  measurements  together  in  obtaining  our 
estimates,  but  the  basic  estimation  process  requires  an  orthgonal  sampling  grid,  which  the 
union  of  the  two  grids  is  not.  Figure  4  illustrates  the  irregular  grid  obtained  by  overlapping 
the  S-N  and  W-E  tracks  of  figures  2  and  3. 

A  direct  approach  to  extending  the  method  so  that  the  survey  data  from  both  files  could 
be  used  is  to  superimpose  the  two  grids  and  fill  in  the  missing  points,  i.e.  the  interiors  of 
the  squares  between  the  tracks,  by  interpolation.  The  new,  densified  grid  would  then  have 
samples  separated  by  1  km  in  both  directions.  Unfortunately,  the  sample  count  of  the  new 
grid,  illustrated  by  Figure  5,  would  thus  be  much  larger  than  either  of  the  two  original 
grids. 

A  much  simpler  method  has  been  devised  whereby  the  K-L  coefficients  are  estimated 
on  each  of  the  two  grids  separately  and  then  combined  to  form  a  single  K-L  coefficient 
field  from  which  the  estimated  gravity  vector  is  computed.  The  two  fields  are  rectangular 
and  overlap.  For  the  overlapping  portion,  the  mean  of  the  two  fields  is  used.  For  non 
overlapping  portions  of  the  two  fields,  whichever  field  is  present  is  used  to  form  the  com¬ 
bined  field.  Otherwise,  zero  is  used.  This  is  called  the  union  method  because  all  of  the  two 
coefficient  fields  were  used,  i.e.  the  (square)  convex  hull  of  their  union.  A  variation  of  this 
method  was  also  employed,  called  the  intersection  method  because  only  the  overlap,  i.e. 
the  (square)  intersection  of  the  two  coefficient  fields  was  used  to  reconstruct  the  gravity 
vector  field  estimates.  The  K-L  coefficient  field  grids  are  illustrated  in  Figures  6  and  7. 
The  two  combination  methods,  the  union  and  intersection ,  are  illustrated  in  Figures  8 
and  9. 
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Figure  5:  Densified  Dual  Grid 


Figure  6:  K-L  Coefficient  Grid:  File  1,  S-N  Tracks 
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Figure  7:  K-L  Coefficient  Grid:  File  2,  W-E  Tracks 


Figure  8:  K-L  Coefficient  Grid:  Combined  Files  1  &  2,  Union  Method 


Figure  D:  K-L  Coefficient  Grid:  Combined  Files  1  &  2,  Intersection  Method 
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2.6  Error  Analysis 

The  primary  error  analysis  tool  employed  is  the  three  dimensional  plot  of  gravity  vector 
estimates  and  truth  values,  which  can  be  compared  visually.  A  error  analysis  based  on 
the  plots  was  undertaken,  in  which  the  minima  and  maxima  of  the  estimate  and  truth 
fields  were  employed  to  compute  some  simple  indicators.  Both  the  plots  and  the  indicators 
appear  in  this  report. 

The  minimum  and  maximum  values  over  the  plotted  region  of  the  estimates  Tx,  Ty,  T:, 
and  the  truth  values  and  Tz  were  taken  from  the  graphs  and  fed  into  an  awk  [1] 

program.  The  following  indicators  were  computed  and  tabulated  for  each  estimate: 

(1) 

(-) 

(3) 


. .  max  Ti  —  max  T*  +  min  7} 
bias,-  = - - - 


•  mill . 


where  i  =  x,  y,  z. 


relativebias,-  = 

maxi,-  —  mini,- 
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3  SOFTWARE  VALIDATION 

3.1  Dipole  Simulation 

The  estimator  software  was  tested  with  deterministic  controlled  data  generated  by  a  simu¬ 
lation  of  the  gravity  field  due  to  a  single  gravity  dipole.  Whether  or  not  such  dipoles  exist 
in  nature  is,  of  course,  entirely  immaterial;  they  were  selected  because  of  their  theoretical 
simplicity,  the  ease  of  computing  them,  and  the  ease  of  interpreting  the  results.  Formu¬ 
las  for  a  dipole  potential  and  its  first  and  second  gradients  were  derived  from  elementary 
physics  and  differential  calculus.  The  gravity  gradient  field  was  computed  at  altitude  and 
the  gravity  vector  field  was  computed  at  ground  level  for  comparison.  A  104  x  104  sample 
grid  with  1  km  was  used.  The  choice  of  104  was  dictated  by  the  fact  that  1)  it  is  even,  which 
makes  the  cosine  matrices  Equations  (97)  and  (99)  invertible  and  their  inverses  expressible 
in  terms  of  the  fast  cosine  transform,  and  2)  104  +  1=105=3  *  5  •  7  is  highly  composite,  a 
property  which  makes  the  fast  sine  and  cosine  transforms  computationally  more  efficient. 
Row  spacing  is  1  km  in  both  directions.  The  dipole  is  horizontal,  centered  in  the  region, 
and  oriented  antiparallel  to  the  x  =  y  line,  i.e.  parallel  to  the  vector  (x,  y,  z )  =  (—1,  —1, 0), 
at  a  depth  of  105/8  km.  The  gravity  gradient  was  computed  at  a  height  of  2  km,  and  the 
gravity  vector  was  computed  at  a  height  at  0  km. 

3.2  Dipole  Math  Model 

The  gravitational  potential  T  due  to  a  point  mass  m  is  given  by 

T  =  mr"1  (4) 

where  r  is  the  distance  from  the  mass  point  to  the  observation  point.  Suppose  we  have  two 
point  masses  +m  and  —  m  at  distances  of  r+  and  r_,  respectively,  from  the  observation 
point.  Then  the  potential  at  the  observation  point  is 

T  =  m(r;1-rl1)  (5) 

If  the  masses  are  separated  by  the  vector  6,  then  linearizing  the  inverse  distance,  this 
becomes 

T  «  m6  •  -^r-1  (6) 

dr 

where  r  «  r+  «  r_  is  the  norm  of  the  vector  r  from  the  midpoint  between  the  two  masses 
to  the  observation  point.  Defining  the  mass  dipole 

p  =  m6  (7) 

and  letting  the  distance  ||<5||  — ►  0  while  preserving  the  product  in  equation  (7),  equation 
(6)  becomes  exact: 

T  =  (8) 
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To  synthesize  the  dipole  field,  we  need  equations  for  the  potential  T ,  the  gravity  vec- 
tor  VT  and  the  gravity  gradient  tensor  VVT.  Switching  to  component  notation,  direct 
differientation  yields: 

T=J2fikdkr~l  (9) 

k 

T,  =  dlT  =  '£(ikdktr-1  (10) 

k 

Tim  =  dimT  =  YLVkdklmr'1  (11) 

k 

Note  that  dki...  is  shorthand  notation  for  the  partied  derivative  operator  with  respect  to 

the  vector  components  rk,  ri, - Also,  all  indices  range  over  the  set  {#,  y,  z}.  The  partial 

differientation  of  powers  or  r  is  facilitated  by  the  formula 

dkrn  =  nrn~2rk 

The  partial  derivatives  are 

di tr-1  =  -r~3rk  (13) 

dkiT~l  =  3r-5rfcr;  —  r~36ki  (14) 

dkimr -1  =  -15r~7r*r/rm  +  3  r_5(r*$fm  +  r/<5mfc  +  rmSki)  (15) 

where  6ki  is  the  Kronecker  delta.  In  the  dipole  simulation  program,  where  correctness  takes 
precedence  over  speed,  the  partials  given  by  equations  (13)— (15)  are  computed  and  sub¬ 
stituted  into  equations  (9)— (11)  yielding  the  gravity  potential,  gravity  vector,  and  gravity 
gradient  tensor. 

For  completeness  of  presentation,  we  note  that  the  substitution  can  be  performed  sym¬ 
bolically  and  simplified  to  obtain  the  following  formulas  for  the  potential  and  its  derivatives: 

T  =  -r^Yh^rk  (16) 

k 

Ti  =  r"3[3r"2(£3  -  /q]  (17) 

k 

Tim  =  3r-5[(]T)  Hkrk)(-5r~2rtrm  +  Sim)  +  rmm  +  rm/q]  (IS) 

k 

These  equations  are  slightly  more  general,  as  well  as  more  compact,  than  equations  (2.1- 
4)-(2.1-13)  of  White  [4],  which  are  restricted  to  vertically  oriented  dipoles.  By  setting 
Hx  =  Hy  =  0  in  equations  (16)— (18),  the  White  equations  may  be  derived,  which  were 
supposedly  used  to  produce  the  NSWC  data  described  in  this  report. 

3.3  Simulation  Results 

3.3.1  Gradient  Measurements 

The  simulated  gravity  gradient  measurements  are  plotted  in  Figures  13-18. 
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3.3.2  Gravity  Truth 

The  simulated  true  gravity  vector  components  are  plotted  in  Figures  10-12. 

3.3.3  Single  Channel 

The  estimated  gravity  vector  components,  based  on  measured  single  components  of  the 
gravity  gradient,  are  plotted  in  Figures  19-42.  Also  plotted  we  the  K-L  coefficient  fields 
a. 


3.3.4  Full  Tensor 

The  estimated  gravity  vector  components,  based  on  the  measured  full  gravity  gradient 
tensor,  are  plotted  in  Figures  43-46.  Also  plotted  is  the  K-L  coefficient  field  a. 

3.4  Error  Analysis 

Visual  comparison  of  the  plots  shows  excellent  agreement  between  the  truth  values  of  the 
gravity  vector  and  the  estimates.  Edge  effects  can  be  seen  clearly  in  the  estimate  plots. 
Table  1  shows  the  results  of  a  simple  error  analysis  based  on  just  the  maxima  and  minima 
of  the  true  and  estimated  fields.  In  this  and  subsequent  tables,  the  first  three  rows  in 
the  table  represent  the  truth  values  themselves,  hence  the  unity  gain  and  zero  bias.  In 
subsequent  tables,  missing  row  entries  correspond  to  estimates  which  exhibited  directional 
instability  (discussed  later  in  the  report)  so  that  the  minimum  and  maximum  values  are 
not  meaningful. 
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Table  1:  Single  Dipole  Gain  and  Bias. 


min 

max 

-9.37e-14 

3.11e-13 

-9.37e-14 

3.11e-13 

-3.79e-13 

3.79e-13 

-9.76e-14 

-1.05e-13 

-3.66e-13 

3.07e-13 

2.92e-13 

3.66e-13 

-1.05e-13 

-9.76e-14 

-3.66e-13 

2.92e-13 

3.07e-13 

3.66e-13 

-9.73e-14 

-9.73e-14 

-3.79e-13 

3.08e-13 

3.08e-13 

3.79e-13 

-9.73e-14 

-9.73e-14 

-3.84e-13 

3.06e-13 

3.06e-13 

3.84e-13 

-9.7e-K 

-1.08e-13 

-3.86e-13 

3.07e-13 

3.09e-13 

3.86e-13 

-1.08e-13 

-9.7e-14 

-3.86e-13 

3.09e-13 

3.07e-13 

3.86e-13 

-9.69e-14 

-9.69e-14 

-3.78e-13 

3.08e-13 

3.08e-13 

3.78e-13 

bias 


0 


0 


0 


4.31e-15 

1.49e-14 

0 


1.49e-14 

4.31e-15 

0 


3.65e-15 

3.65e-15 

0 


4.37e-15 

4.38e-15 

-6e-18 


3.65e-15 

8.21e-15 

-8e-18 


8.22e-15 

3.64e-15 

-6e-18 


3.18e-15 

3.18e-15 

-2.01e-18 


relbias 


0 


0 


0 


0. 

0. 

0 


0. 
0. 
0 


0.00902 

0.00902 

0 


-7.91e-06 


-1.05e-05 


0.0203 

0.00899 

-7.91e-06 


0.00785 

0.00785 

-2.G4e-06 


4  NSWC  DATA 

4.1  Data  Tape 

The  algorithm  was  tested  on  foreign  data  obtained  from  the  Naval  Surface  Weapons  Center 
(NSWC)  [3],  The  data  were  synthesized  from  a  multilayer  vertical  dipole  array.  The 
magnitude  of  the  dipoles  was  randomly  selected  in  such  a  way  as  to  produce  gravity  fields 
which  are  statistically  similar  to  those  observed  in  a  certain  region  of  Texas.  Details  of  the 
model  are  found  in  White  [4], 

The  simulation  data  were  supplied  on  a  magnetic  tape  which  we  refer  to  as  “tape  l‘\ 
There  are  three  files  of  data  on  this  tape:  the  first  two  contain  simulated  gravity  gradient 
measurements  such  as  would  be  obtained  by  a  survey  aircraft.  The  third  contains  the 
gravity  force  vector  at  zero  altitude.  File  1  data  are  presented  as  south-to-north  (S-N) 
tracks,  file  2  date  are  presented  as  west-to-east  (W-E)  tracks,  and  file  3  data  are  presented 
at  the  intersections  of  the  S-N  and  W-E  tracks.  The  sampling  density  is  greater  along 
track  than  across  track. 

4.2  Grid — Index  Mapping 

Each  file  was  divided  into  records,  each  record  constituting  a  sample  point.  The  x  (east) 
and  y  (north)  coordinates  for  the  point  were  given  and  then  either  the  six  components  of  the 
gravity  gradient  tensor  (files  1  and  2)  or  the  three  components  of  the  gravity  disturbance 
vector  (file  3). 

The  range  of  x  and  y  was  —250  km  to  +250  km.  In  the  direction  of  the  tracks,  the 
sampling  interval  was  1  km  and  in  the  perpendicular  direction  it  was  5  km.  This  makes  the 
full  arrays  501  x  101  points  (file  1),  101  X  501  points  (file  2)  and  101  x  101  points  (file  3). 

Each  tape  file  was  read  into  the  computer  and  processed.  Each  physical  quantity 
represented  was  placed  in  a  matrix  whose  elements  represented  the  sampling  nodes  in  the 
grid.  The  value  of  the  coordinate  x  or  y  was  mapped  linearly  as  appropriate  so  that  the 
lowest  index  value  of  1  represented  the  lowest  coordinate  value,  and  the  highest  index  value 
represented  the  highest  coordinate  value.  The  matrices  were  then  output  sequentially;  this 
permitted  the  estimation  software  to  be  structured  so  that  one  measurement  field  at  a 
time  could  be  processed.  This  sequential  processing  of  the  measurements  greatly  reduces 
the  amount  of  storage  needed  to  perform  estimation  from  multi-component  (vector  and 
tensor)  measurements. 

4.3  Cropping  of  the  Data  Region 

The  number  of  samples  in  each  direction  recorded  on  the  tape  was  selected  without  regard 
to  the  requirements  of  any  particular  estimation  algorithm.  Because  the  cosine  matrix 
equations  (58)-(59)  must  have  an  even  number  of  rows  and  columns  in  order  to  be  invert¬ 
ible,  we  require  an  even  number  of  samples.  Furthermore,  in  order  to  use  fast  sine  and 
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cosine  transforms  efficiently,  the  number  of  samples  in  each  direction  should  be  one  less 
than  a  highly  composite1  (odd)  number.  Of  the  101  available  points  in  the  cross-track 
direction,  we  therefore  used  98  =  32  •  11  —  1  corresponding  to  the  range  -240,  —235,  . . . , 
+245  km.  Of  the  501  available  points  in  the  along  track  direction,  we  used  494  =  32>5T1  — 1 
corresponding  to  the  range  -246,-245,  . . .  ,+247  km.  Table  2  summarizes  the  grid  infor¬ 
mation.  The  symbols  Ax,  and  Ay  denote  the  sampling  intervals  along  each  axis  (east  and 
north),  K  and  L  denote  the  numbers  of  samples,  and  A  =  (K  +  l)Ax  and  B  =  (L  +  l)Aj/ 
denote  the  region  dimensions. 


4.4  Data  Preprocessing — Transformation 

The  data  were  presented  in  north-east-down  (NED)  coordinates.  They  were  originally 
believed  to  have  been  in  east-north-up  (ENU)  coordinates.  Because  our  analysis  and 
software  were  based  on  ENU  (really  any  right-handed  frame  with  the  z  axis  pointing 
upward)  we  applied  a  transformation  from  NED  to  ENU  to  all  the  data.  The  equations 
for  this  transformation  are: 


/  J 

®  i  V  >  2 

II 

~l 

M 

(19) 

rpt  rpt  rpt 
**  y )  xz 

—  Ty ,  ,  Tz 

(20) 

/TV  rpt  rpt 

*  *  yy)  ±  zz 

=  Tyy  )  Txx  ,  Tzz 

(21) 

rpt  rpt  rpt 
±xy)  xyz)  xxz 

=  TXy)  -~TXzy  —TyZ 

(22) 

where  unprimed  coordinates  are  in  the  old  NED  frame  and  the  primed  coordinates  are  in 
the  new  ENU  frame.  Here  also  T  is  the  disturbing  potential  and  the  x,  y,  and  z  subscripts 
denote  partial  differientation. 

Incidentally,  the  transformation  works  equally  well  in  reverse. 


'The  product  of  small  primes. 


Table  2:  Sampling  Grids 


file  no. 

Ax  (km) 

Ay  (km) 

A  (km) 

B  (km) 

K 

a 

1 

1 

5 

495 

495 

494 

98 

2 

5 

1 

495 

98 

494 

3 

5 

5 

495 

495 

98 

98 
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4.5  Gradient  Measurements 

4.5.1  File  1 


The  simulated  gravity  gradient  measurements  along  S-N  tracks,  recorded  on  tape  1,  file  1, 
are  plotted  in  Figures  47-52. 

4.5.2  File  2 

The  simulated  gravity  gradient  measurements  along  W-E  tracks,  recorded  on  tape  1,  file  2, 
are  plotted  in  Figures  81-86. 

4.6  Gravity  Truth,  File  3 

The  simulated  true  gravity  vector  components,  recorded  on  tape  1,  file  3,  are  plotted  in 
Figures  123-125. 
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5  TEST  RESULTS— NSWC  DATA 

5.1  Single  Channel,  File  1 

The  estimated  gravity  vector  components,  based  on  measured  single  components  of  the 
gravity  gradient  (tape  1,  file  1)  measured  on  S-N  tracks,  are  plotted  in  Figures  53-76.  Also 
plotted  are  the  K-L  coefficient  fields  at.  Table  3  shows  the  gain  and  bias  indicators  for  this 
scries  of  estimates. 


5.2  Full  Tensor,  File-  1 

The  estimated  gravity  vector  components,  based  on  the  measured  full  gravity  gradient 
tensor  (tape  1,  file  1)  measured  on  S-N  tracks,  are  plotted  in  Figures  77-80.  Also  plotted 
is  the  K-L  coefficient  field  d.  Table  3  shows  the  gain  and  bias  indicators  for  this  series  of 
estimates. 


Table  3:  File  1  Gain  and  Bias. 


meas 

est 

min 

max 

gain 

bias 

relbias 

Tx 

Tx 

-38.2 

45 

1 

0 

0 

Ty 

Ty 

-24.5 

59.6 

1 

0 

0 

Tz 

Tz 

-53.4 

64.3 

1 

0 

0 

^xx 

Tx 

-45.1 

52.4 

1.17 

-0.229 

-0.00276 

Tyy 

Tx 

-37 

32.3 

0.833 

5.76 

0.0692 

Tyy 

Ty 

-44.3 

35.9 

0.953 

21.7 

0.258 

y 

Tz 

-66.4 

53.7 

1.02 

11.8 

0.1 

Tgz 

Tx 

-41.4 

40.9 

0.989 

3.69 

0.0444 

T„ 

Ty 

-47.3 

34.9 

0.976 

23.8 

0.283 

Tsz 

Tz 

-59.8 

59.6 

1.02 

5.56 

0.0473 

TXy 

Tx 

-40.5 

40.2 

0.971 

3.57 

0.043 

^yz 

Tx 

-41.8 

40.8 

0.993 

3.92 

0.0472 

T 

±yz 

Ty 

-47.3 

35.5 

0.9S4 

23.4 

0.279 

71 

1yz 

Tz 

-77.1 

65.5 

1.21 

11.3 

0.0959 

Txs 

Tx 

CSJ 

i — i 

* 

46.5 

1.05 

0.776 

0.00933 

full 

Tx 

-40.3 

41.4 

0.983 

2.89 

6.0348 

full 

Ty 

-47 

31.9 

0  938 

25.1 

0.298 

full 

Tz 

-56.9 

59.4 

0.989 

4.16 

0.0354 
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5.3 


Single  Channel,  File  2 

The  estimated  gravity  vector  components,  based  on  measured  single  components  of  the 
gravity  gradient  (tape  1,  file  2)  measured  on  W-E  tracks,  are  plotted  ir  Figures  87-110. 
Also  plotted  are  the  K-L  coefficient  fields  a.  Table  4  shows  the  gain  and  bias  indicators 
for  this  series  of  estimates. 

5.4  Full  Tensor,  File  2 

The  estimated  gravity  vector  components,  based  on  the  measured  full  gravity  gradient 
tensor  (tape  1,  file  2)  measured  on  W-E  tracks,  are  plotted  in  Figures  111-114.  Also 
plotted  is  the  K-L  coefficient  field  a.  Table  4  shows  the  gain  and  bias  indicators  for  this 
series  of  estimates. 


Table  4:  File  2  Gain  and  Bias. 


rneas 

est 

min 

max 

gain 

bias 

relbias 

Tx 

Tx 

-38.2 

45 

1 

0 

0 

Ty 

Ty 

-24.5 

59.6 

1 

0 

0 

Tz 

Tz 

-53.4 

64.3 

1 

0 

0 

TXx 

Tx 

O 

i 

42.5 

0.991 

2.18 

0.0262 

Txx 

Ty 

-4'  .7 

29.8 

0.861 

24 

0.286 

'I'xx 

Tz 

-52.5 

62.4 

0.977 

0.495 

0.00421 

Tyy 

Ty 

-52.4 

34.8 

1.04 

26.4 

0.314 

Tzz 

Tx 

-39.3 

39.6 

0.949 

3.26 

0.0392 

Tzz 

Ty 

-46.6 

32.8 

0.943 

24.4 

0.291 

Tzz 

Tz 

-58.7 

57.9 

0.991 

5.85 

0.0497 

Txy 

Ty 

-50 

37.9 

1.04 

23.7 

0.281 

Ty . 

Ty 

-47.6 

36.2 

0.996 

23.2 

0.276 

Txz 

Tx 

-40.1 

43,1 

1 

1.91 

0.023 

Txz 

Ty 

-47.6 

40.2 

1.04 

21.3 

0.253 

Txz 

Tz 

-52.8 

72 

1.06 

-4.14 

-0.0352 

full 

Tx 

-39.4 

41.4 

0.972 

2.44 

0.0294 

full 

Ty 

-46 

31.2 

0.917 

24.9 

0.296 

full 

Tz 

-55.4 

59.3 

0.975 

3.54 

0.0301 
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5.5  Full  Tensor,  Combined  Files 

5.5.1  Intersection  Method 

The  estimated  gravity  vector  components  and  K-L  coefficients,  based  on  using  the  inter¬ 
section  method ,  on  the  measured  full  gravity  gradient  tensor  from  both  S-N  tracks  (tape  1, 
file  1)  and  W-E  tracks  (tape  1,  file  2),  are  plotted  in  Figures  115-118.  Table  5  shows  the 
gain  and  bias  indicators  for  this  series  of  estimates. 


Table  5:  Files  1  &  2  (intersection  method)  Gain  and  Bias. 


meas 

est 

min 

max 

gain 

bias 

relbias 

Tx 

Tx 

-38.2 

45 

1 

0 

0 

Ty 

Ty 

-24.5 

59.6 

1 

0 

0 

% 

Tz 

-53.4 

64.3 

1 

0 

0 

full 

Tx 

-39.1 

40.3 

0.955 

2.86 

0.0344 

full 

Ty 

-45.2 

31.2 

0.908 

24.6 

0.292 

full 

Tz 

-55.1 

59.5 

0.974 

3.27 

0.0278 

5.5.2  Union  Method 

The  estimated  gravity  vector  components  and  K-L  coefficients,  based  on  using  the  union 
method,  on  the  measured  full  gravity  gradient  tensor  measured  on  both  S-N  tracks  (tape  1, 
file  1)  and  W-E  tracks  (tape  1,  file  2),  are  plotted  in  Figures  119-122.  Table  6  shows  the 
gain  and  bias  indicators  for  this  series  of  estimates. 


Table  6:  Files  1  &  2  (union  method)  Gain  and  Bias. 


meas 

est 

min 

max 

gain 

bias 

relbias 

Tx 

-38.2 

45 

i 

0 

0 

Ty 

wa 

-24.5 

59.6 

i 

0 

0 

Tz 

-53.4 

64.3 

i 

0 

0 

EH 

Tx 

-40.5 

1 

m 

3.13 

i 

0.0376 

m 

Ty 

-46.4 

0.921 

25.2 

0.299 

full 

ULj 

|  -57.6 

60 

1 

4.26 

0.0362 
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ANALYSIS  OF  RESULTS 


6.1  Directional  Instability 

An  interesting  phenomenon  which  occurred  with  some  single  channel  measurement  esti¬ 
mates  is  directional  instability,  i.e.  the  appearance  of  a  surface  which  is  smooth  when 
scanned  along  one  horizontal  axis  but  which  oscillates  wildly  when  scanned  in  the  direc¬ 
tion  of  the  other  horizontal  axis.  The  occurrence  of  these  instabilities  has  been  noted  from 
the  graphs  and  tabulated  for  file  1  (S-N  tracks)  and  file  2  (W-E  tracks)  single  channel 
estimates.  We  observe  that  only  the  gravity  vector  component  which  can  be  consistently 
estimated  is  the  horizontal  one  at  right  angles  to  the  tracks,  that  is  in  the  sparsely  sam¬ 
pled  direction.  Only  half  the  remaining  combinations  yield  smooth,  reasonable  estimates. 
The  other  half  (six  combinations  for  each  file)  are  unstable  in  the  direction  parallel  to  the 
tracks,  i.e.  the  direction  on  the  higher  sampling  density.  For  file  1  we  can  state  the  rule 
that  an  estimate  will  be  unstable  in  the  y  (dense)  direction  if  and  only  if  the  measurement 
is  an  x  partial  derivative  and  the  estimated  quantity  is  not.  For  file  2  the  same  rule  applies 
if  only  we  exchange  the  x  and  y  in  the  rule  for  file  1.  Thus  a  gravity  component  in  either 
the  densely  sampled  or  vertical  direction  can  not  be  obtained  from  a  gradient  component 
in  the  sparsely  sampled  direction. 

More  investigation  is  required  to  determine  whether  the  problem  of  directional  instabil¬ 
ity  for  single  channel  measurements  with  unequal  sampling  intervals  in  the  two  horizontal 
directions,  is  fundamental  or  can  be  overcome.  In  any  case  they  do  not  pose  a  serious 
problem  so  long  as  full  tensor  measurements,  which  do  no  exhibit  directional  instabiliity , 
are  available.  Even  if  the  full  tensor  is  not  available,  the  method  is  general  enough  to 
accomodate  any  subset  of  the  six  (five  independent)  gradient  tensor  components,  which 
presumably  would  be  sufficient  to  overcome  the  instability. 


6.2  Gain  and  Bias 

In  examining  the  tabulated  estimator  gain  and  bias  errors,  we  observe  that  single  axis  gains 
ranged  from  a  decrease  of  17%  to  an  increase  of  21%  for  file  1  and  from  a  decrease  of  14% 
to  an  increase  of  6%  for  file  1.  Relative  biases  vary  from  —3%  to  +30%  for  file  1  and  from 
—4%  to  +31%  for  file  2.  Gain  and  bias  errors  based  on  full  tensor  gradient  measurements 
show  much  less  spread,  however,  than  the  single  axis  errors:  the  gains  range  from  6%  less 
than  unity  (file  2)  to  1%  less  than  unity  (file  1).  Relative  biases  range  from  +3%  to  +30% 
(both  files).  Interestingly,  estimates  of  Ty  show  a  consistent  relative  bias  of  approximately 
30%.  Results  are  similar  for  the  combined  estimates  using  both  the  intersection  and  union 
methods,  although  the  latter  seems  slightly  better  on  gain  and  the  former  slightly  better 
on  relative  bias. 

The  slight  attenuation  typical  of  most  of  the  estimates  is  fairly  easily  explained  by 
noting  that  the  estimator  is  a  low  pass  filter.  This  is  because  the  estimator  knows  that  the 
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measurement  noise  is  a  wide  band  signal.  By  filtering  out  the  higher  spatial  frequencies, 
the  estimator  will  dull  the  sharp  peaks  which  more  than  likely  include  the  minima  and 
maxima,  resulting  in  a  smaller  peak-to-peak  magnitude  and  hence  a  less  than  unity  gain 
by  our  definition. 


6.3  Boundary  and  Edge  Effects 

Likely  error  sources  are  b  j  ndary  and  edge  effects.  By  the  former  we  mean  errors  resulting 
from  assuming  boundary  conditions  which  do  not  in  fact  obtain.  Boundary  effects  which 
are  in  evidence  only  near  the  edges  of  the  region  are  known  as  edge  effects.  The  potential 
model  of  a  pure  sine  Fourier  series  in  x  and  y  given  by 

00  CO 

T{x ,  y,  z)  -  Yj  Yj  sin(ajk®)  sin (fyy)  exp(-cktz)  (23) 

Jt=i  j=i 

where 

ak  =  kir/A,  bi  =  Iw/B,  ck\  =  \Ja\  +  bj  (24) 

implies  certain  boundary  conditions:  1)  T  is  periodic  with  period  (2Ay2B)>  i.e. 

T(x,  y,z )  =  T(  x  +  2A,  y,  z)  =  T(.r,  y  +  2  B,  z)  (25) 

and  2)  T  is  odd  symmetric  in  x  and  y,  i.e. 

T(.r,  y,  z)  =  —T(—x,  y,  z)  =  ~T(x,  -y,  z)  (26) 

These  follow  directly  from  the  properties  of  the  sine  function  and  are  readily  verifiable  from 
the  representation  of  T(x,y,z)  given  by  equation  (23).  A  corollary  to  these  conditions  is 
that  T  is  odd  symmetric  about  any  of  the  survey  region  boundaries: 

T(A  +  y,  z)  s  -T(A  -  x,  y,  z),  T(x,  B  +  y,  z)  =  -T(x,  B  -  y,  z)  (27) 

It  follows  that  the  mass  distribution  which  induces  T  must  satisfy  the  same  conditions  as 
the  potential,  i.e.  each  mass  particle  or  feature  will  have  a  two-dimensional  infinite  series 
of  reflective  and  periodic  images  of  itself,  some  of  which  are  sign  reversed. 

To  the  extent  that  the  real  mass  distribution  meets  these  conditions,  the  estimator 
will  produce  correct  results.  Generally,  however,  the  mass  distribution  will  not  meet  these 
symmetry  and  periodicity  conditions.  In  such  cases,  good  estimates  may  still  be  possible 
away  from  the  boundaries.  The  sources  of  error  may  be  masses  outside  the  region,  which 
are  necessarily  not,  modeled,  and  the  reflective  and  periodic  images  located  outside  the 
region  which  are  unavoidably  and  erroneously  modeled.  Since  the  potential  due  to  a.  mass 
particle  declines  with  distance  from  the  observer,  estimates  can  be  expected  to  be  good 
near  the  surface  and  away  from  the  edges,  provided  there  are  not  unduly  significant  mass 
distributions  near  or  beyond  the  edges  or  at  great  depth  compared  to  A  and  B. 
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Let  us  remark  that  the  symmetry  condition  can  be  done  away  with  by  switching  to  a 
potential  representation  which  includes  the  cosine-cosine  as  well  a  mixed  sine  and  cosine 
terms  of  the  Fourier  series 

T(*,S,*)  =  t  (28) 

CO  CO 

X  X  aVi  sin(afc»)  sin(6,j/) 

/=i 

00  CO 

+  sin(ajfca:)  cos(6/t/) 

Jt=i  l=o 

CO  CO 

+  X  X)  aJW  cos(akx)  sm(b,y) 
k=o  /=i 

CO  oo 

+  X  X  aVt  cos(ajtx)  cos(bty) 

k= 0  1=0 

]  exp(— C«2 ) 

or  the  equivalent  complex  Fourier  series.  This  alternative  approach  requires  data  to  be 
expanded  to  a  full  period,  i.e.  a  rectangle  of  dimensions  2A  x  2 B.  It  has  the  advantage  of 
being  able  to  handle  non-zero  potentials  on  the  boundary. 

An  alternative  approach  to  dealing  with  non-zero  boundary  potentials  is  to  fit  a  simple 
harmonic  field  to  the  region  boundaries  and  then  apply  the  existing  estimation  algorithm  to 
estimating  the  residual  field,  which  will  necessarily  meet  the  zero  boundary  value  condition. 
The  final  field  estimate  is  the  sum  of  the  simple  field  which  fits  the  boundaries  and  the 
estimated  residual  field.  Such  a  fit  must  be  done  using  measurements  of  the  gravity  gradient 
on  the  region  boundary.  Currently  these  measurements  are  not  utilized  in  our  algorithm. 
A  model  for  the  boundary-fitting  harmonic  field  is  obtained  from  the  full  Fourier  series 
equation  (28)  by  restricting  the  range  of  the  cosine  summation  indices  to  a  maximum  of  1. 

T(x,ytz)  =  [  (29) 

X  X  <*ki  sin  (a**)  sin  (%) 

tel  /=l 
co  1 

+  EXqh  sin(a*s)  cos  (bty) 

k=l  1=0 
1  co 

+  x  X  ati  c°s(akx)  sin(biy) 

k= 0  /=1 
1  1 

+  XEaS  cos  (akx)  cos(biy) 
k=o  1=0 

]exp  (-cmz) 

The  details  of  estimating  the  boundary  coefficients  are  beyond  the  scope  of  the  present 
research.  The  purpose  of  introducing  them  here  is  to  demonstrate  a  model  for  the  errors 
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incurred  when  the  boundary  values  of  the  potential,  being  observed  through  the  gradients, 
are  non-zero. 

The  Fourier  coefficients  in  equation  (29)  tare  obtained  sequentially.  The  a£j  ar  deter¬ 
mined  first  so  that  the  four  cosine-cosine  terms  produce  a  harmonic  field  which  tits  the 
corners  of  the  potential.  The  residual  field  is  computed  by  subtracting  the  cosine-cosine 
terms  from  the  total  field  T(x,  y,  z)  leaving  be  a  harmonic  field  with  zero  potential  on  the 
four  corners.  Then  each  set  of  mixed  coefficients  is  computed  to  fit  the  interiors  of  one 
parallel  set  of  sides.  A  new  residual  field  is  computed  which  has  homogenous  boundary 
conditions  and  is  thus  representable  in  terms  of  sine-sine  terms  only.  But  we  have  already 
solved  the  problem  of  estimating  the  sine-sine  coefficients  by  the  Karhunen-Loeve  method. 
The  details  of  this  multistep  estimation  process  based  on  equation  (29)  are  beyond  the 
scope  of  this  research.  Note,  however,  that  besides  providing  a  means  of  extending  the 
K-L  method  to  handle  non-homogenous  problems,  it  gives  one  a  model  for  determining 
what  types  of  errors  may  be  present  when  the  homogenous  K-L  method  is  applied  to 
non-homogenous  data  as  was  done  in  this  study. 

By  visualizing  several  of  these  boundary  harmonics,  the  reader  can  imagine  what  types 
of  errors  may  be  present.  Some  examples: 

•  constant  in  both  x  and  y  directions 

•  cosine  in  *,  constant  in  y 

•  cosine  in  a:  and  y 

•  arbitrary  Ci  curve  in  a:  vanishing  at  the  endpoints,  constant  in  y 

•  arbitrary  Ci  curve  in  y  vanishing  at  the  endpoints,  constant  in  x 

By  this  analysis  we  demonstrate  a  rather  large  family  of  harmonic  surfaces  from  which  to 
select  one  which  corrects  for  the  boundary  values. 
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SUMMARY,  CONCLUSIONS,  RECOMMENDA¬ 
TIONS 


7.1  Summary 

The  Karhunen-Loeve  algorithm  was  implemented,  validated  on  controlled  data,  and  tested 
on  foreign  data.  Successful  recovery  of  the  unknown  gravity  field  was  accomplished  in 
most  cases,  the  only  exceptions  being  certain  single-component  measurement  cases  where 
the  sampling  interval  is  different  in  the  two  horizontal  directions.  Plots  of  the  measured 
gravity  gradients  and  of  the  true  and  estimated  gravity  vectors  have  been  presented.  A 
capability  to  simulate  dipole  fields  was  developed  and  applied  to  produce  controlled  data 
based  on  a  single  horizontal  gravity  dipole;  the  single  dipole  data  were  used  for  software 
validation. 

Error  analyses  were  based  on  qualitative  comparison  of  the  plots  and  on  quantitative 
comparison  of  indicators  computed  from  the  minima  and  maxima  of  the  plotted  gravity 
vector  fields.  Since  a  method  designed  for  homogeneous  potential  fields,  i.e.  zero-valued  on 
the  rectangular  survey  region  boundaries,  was  being  applied  to  data  derived  from  a  non- 
homogenous  potential,  the  estimate  regions  were  cropped  on  the  hypothesis  that  boundary 
and  edge  effect  errors  would  be  mostly  restricted  to  near  the  rectangle  boundaries.  This 
hypothesis  was  at  least  partially  justified  by  the  results.  Ways  to  extend  the  method  to 
deal  properly  with  the  non-homogenous  boundary  conditions  were  suggested. 

In  the  process  of  implementing  the  K-L  algorithm  in  software,  several  technical  advances 
were  made.  These  include  1)  an  arrangement  of  equations  which  minimizes  computer  stor¬ 
age  of  arrays  and  permits  (field)  measurements  to  be  made  sequentially;  2)  the  application 
of  fast  sine  and  cosine  transforms  to  speed  computation  dramatically,  especially  for  large 
sample  grids. 


7.2  Conclusions 

We  conclude  that  the  Karhunen-Loeve  method  of  gravity  estimation  from  gravity  gradient 
measurements  as  implemented  for  this  study  is  a  viable  estimation  technique.  The  present 
experiments  did  not  attempt  to  account  for  non-homogenous  boundary  values,  which  lim¬ 
ited  the  region  of  acceptable  estimates  to  center  of  the  measurement  region,  i.e.  away 
from  the  boundaries.  It  is  believed  that  better  performance  over  the  whole  region  will 
be  attained  when  the  method  is  enhanced  in  the  ways  recommended.  Still  the  method 
worked  well  within  its  present  limitations.  Because  of  its  excellent  results,  computational 
economy,  and  suitability  for  small  computers,  its  eventual  wide  adoption  is  suggested. 
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7.3  Recommendations 

The  principal  recommendation  for  further  work  is  that  the  method  be  extended  to  account 
for  non- homogeneous  boundary  conditions.  Two  approaches  have  been  pointed  out:  first, 
adding  a  sufficient  number  of  cosine  and  mixed  sine-cosine  harmonics  to  the  disturbing 
potential  model  and  estimating  the  additional  coefficients  from  the  boundary  data  alone; 
and  alternatively,  by  adding  the  full  complement  of  cosine  and  mixed  sine-cosine  terms 
and  expanding  the  data  region  to  a  full  period  in  the  x  and  y  directions.  A  variation  on 
the  latter  method  would  be  to  use  a  complex  Fourier  series  instead  of  the  real  sine-cosine 
series  to  represent  the  potential.  Additional  tests  with  simulated  measurement  noise  are 
recommended. 

We  note  that  the  algorithm  is  not  limited  to  simply  obtaining  estimates  of  the  gravity 
vector  from  measurements  of  the  gradient.  With  minimal  modification  to  the  software, 
combinations  of  measurements  of  the  potential  derivates  of  any  or  different  orders  can  be 
processed  and  likewise  derivatives  of  any  order  can  be  estimated. 

The  mathematical  model  for  the  method  tested  here  can  be  extended  in  potentially 
useful  ways.  For  example,  instead  of  representing  the  source  of  the  disturbing  potential  by 
a  single  white  noise  layer,  multiple  layers  should  be  considered.  These  would  be  placed  at 
various  depths  as  in  the  model  used  to  produce  the  NSWC  data. 
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A  Algorithm  Summary 

A.l  Local  Region 

A  local  region  of  interest  is  represented  as  such,  where  x  and  y  coordinates  are  on  the 
surface  of  the  earth  and  the  z  coordinate  is  vertical  above  the  surface  of  the  earth.  The 
local  region  of  interest  is  specified  by  the  spatial  domain  D  given  as 

j D(x,y)  =  {x,y;  0  <  x  <  A,  0  <  y  <  B]  (30) 

where 


A  -  length  of  the  survey  region 
B  -  width  of  the  survey  region 

The  survey  measurements  are  taken  at  an  altitude  H  above  the  surface  of  the  earth. 


A. 2  Sensor  Signal  Model 

The  gradiometer  sensor  signal  S(x ,  y,  z )  has  six  (6)  components:  three  (3)  inline  signals 
and  three  crossline  signals.  The  six  (6)  signals  of  the  gradiometer  are  given  below 

d2T 


Si(x,y,z) 
S2(x,y,z) 
S3(x,y,z) 
Sa( x,  y,  z) 
S5(x,y,z) 
S6(x,y,z) 


dx2 
!PT 
dy2 
d2T 
dz 2 
cPT 
dxdy 
d*T 
dydz 

cPT 

dzdx 


(31) 

(32) 

(33) 

(34) 

(35) 

(36) 


A. 3  Measurement  Model 

The  measurement  Z(x,y)  at  z  =  H  of  the  gradiometer  signal  S(x,y )  is  contaminated  with 
an  additive  noise  V(x,y)  such  that 

Z(X,y)  =  S(x,y)  +  V(Xiy)  (37) 

The  autocorrelation  function  of  the  additive  noise  is  assumed  to  be  a  diagonal  matrix  such 
that 

Rvv  (xuyi;x2,y2)  =  {V(xuyl)VT(x2, y2)}  (38) 
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a\  0  0  0  0  0 

0  cj\  0  0  0  0 

0  0  a\  0  0  0  £t  sc,  ,  ,on. 

0  0  0  <7?  0  0  «(*>-**)%' -*)  <30> 

0  0  0  0  o\  0 

0  0  0  0  0  4 

The  vector  signal  S(x,y)  and  measurement  noise  V(x,y)  are  assumed  to  be  uncorrelated, 
i.e. 

E  [S(xuyi)VT(x2,y2)}  =  {V(a;i ,  y\ )ST(x2,y2)}  =  0  (40) 

A. 4  Measurement  Grid 

In  order  to  obtain  the  vector  signal  estimates  for  all  the  spatial  points  of  interest  a  definition 
of  the  two-dimensional  grid  is  necessary.  For  the  present  development  the  two-dimensional 
grid  of  data  points  is  defined  to  be  equally  spaced  in  each  direction  with  K  data  points  in 
the  x  -  direction  and  L  data  points  in  the  y  -  direction.  With  Ax  and  Ay  the  grid  spacing 
in  the  x  and  y  directions  respectively  the  spatial  domain  D  is  given  as 


(K  +  l)Ax  =  A 

l<k<I( 

(41) 

(L  +  l)Ay  =  B 

1  <1<L 

(42) 

A. 5  Measurement  Matrices 

'  Zi(xuyi)  Z\(x\ ,  1J2)  •••  Zi(xuyL) 

[Zi]  _  Zi(x2,yi)  : 

KxL  ;  : 

.  Z\{xk -,yi)  . Zi(xk,Vl) 

Z2(xuyi)  Z2(xx,y2)  •••  Z2(xx,yL) 
[Z2]  _  Z2(x2,yi)  : 

KxL  ~  :  : 

.  Z2(xk,Vi)  . Z2(xK,yL) 

'  H*un)  Zz(x\ , y2)  •••  Z3(x, , yL) 

[Z3]  _  Zs(x2,yi)  : 

KxL  :  : 

_Z-j(x  Ktyi)  . Za{xKtyit) 


(43) 


(44) 


(45) 
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(46) 


Z4(X\,yi) 

IT 

to 

*  Z4(xi,yi) 

fa]  = 

Z4{x2,yi) 

| 

KxL 

.  Z4(xKtyi) 

.  .  * 

•  Z\{XK,yL)  . 

'  Zs(xi,yi) 

Zs(xuy2)  • 

•  Z6(xi,yL)  ' 

[Zs]  m 

Zs(x2,yi) 

• 

KxL 

.  Z5(xK,yi) 

••• 

•*  25(®K,  Vl)  . 

Z&(x\,y\) 

Z6(xt ,  2/2)  • 

Z6(xuyL)  ' 

[Ze]  = 

Z6(x2,yi) 

J 

KxL 

.  Z6(xK,yi) 

...  .  , 

••  Z6(xk,Vl)  . 

(47) 


(48) 


A. 6  Transformed  Measurements 


[Zi] 

AB 

[SAA]r[Fi][SAF] 

(49) 

KxL 

(K  +  1)(L  + 1) 

KxK  KxL  LxL 

&) 

AB 

[SAArHZ2][SAr] 

(50) 

KxL 

(K  +  l)(L  +  l) 

KxK  KxL  LxL 

1Z3] 

AB 

[S'AA]t[Z3][5AF] 

(51) 

KxL 

(K  +  1)(L  +  1) 

KxK  KxL  LxL 

[z4] 

AB 

[I  +  2C'][C'AA]T[Z4][C,AF][I  +  2  C) 

(52) 

KxL 

(K  +  1  )(L  + 1) 

KxK  KxK  KxL  LxL  LxL 

[Zs] 

AB 

[SAX]T[Z5}[CAY][I  +  2C] 

(53) 

KxL 

(K  +  1)(L  +  1) 

KxK  KxL  LxL  LxL 

[Zs] 

AB 

[I  +  2C}[CAX]T[Z6}[SAY] 

(54) 

KxL 

(K  +  1  )(L  +  1) 

KxK  KxK  KxL  LxL 

.  L  are  even,  [I]  is  the  K  x  K  or  L  x  L  identity  matrix,  and  [C]  is  a 

K  x  K  or 

LxL  Toeplitz  Circulant  matrix  of  alternating  one’s  (1)  and  zeroes  (0),  a  4  X  4  example 
of  which  is 

"  1  0  10 
0  10  1 
10  10 
0  10  1 


[C}  = 


(55) 
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The  transformation  matrices  are  given  by 


[SAX(iJ)]  =  sin  n(i  x  j)~j~ 

AY 

A  X 

[CAA^i,;')]  =  cos7r (i  x  j)~j~ 

AY 

[CAF(i,  j)]T  =  cos 7r(t  x  ;)— 


1  <  i  <  K 
1  <  j  <  K 

1  <i<L 
l<j<L 

1  <  *  <  Ii 

l<j<  K 

1  <i<L 
1  <j<L 


(56) 

(57) 
(68) 
(59) 


A. 7  Karhunen-Loeve  Coefficient  Estimates 

[^mn]  _  Pmn i  0  £i(ra,  n)  0  7mnl  +  ^mn2  0  Z2(m,  n)  0  7mn2 
AxL  AxL  AxL  AxL  AxL  AxL  AxL 

+/?mn3  ®  £3(m,  «)  0  7m„3  +  ()mn4  0  Z4(m,  n)  0  7mn4 
AxL  AxL  AxL  AxL  AxL  AxL 

+/Smn5  0  z5(m,  n)  0  7mn5  -f  /?mn6  0  Z6(m,  n)  0  7mn6 
AxL  A’xL  A'xL  AxL  A'xL  AxL 


where  0  is  point-by-point  matrix  multiplication  and  not  row-by-column  matrix  multipli¬ 
cation  and, 


n|=1<r?  +  (Ef=1 

(61) 

and 

\  _ 

Amn  “  TiTo 

°mn 

(62) 

and 

^  w  n*  e-cmn|//+X?| 

7mnl  m 

(63) 

/y  „  —  ?  j2.— cmn|A+D| 

7mn2  v/AB  n 

(64) 

<Y  o  —  ^  />2  g-cmn|A+D| 

7mn3  n/ab  mn 

(65) 

7mn4  =  ^/XBambne  Cmn|//+£)| 

(66) 

T  ,  _  -2sgn(if  +  D)  -CmnlA+DI 

7*nn5  *“  \/]45 

(67) 

T  .  _  -2sgn(B  +  B)  Ctmt|//+fli 

7mn6  —  \J  A  B 

(6S) 
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A. 8  Interpolation  Grid 

The  signal  estimates  are  available  at  any  altitude  h.  The  estimates  of  the  Karhunen-Loeve 
coefficients  dmn,  are  on  the  measurement  grid  K  x  L.  For  the  signal  estimates,  a  finer 
interpolated  grid  can  be  used  that  is  defined  by 

(P  +  1)EX  =  A  EAr  <  AX  (69) 

(Q  +  1)E Y  =  B  EF  <  AY  (70) 

A. 9  Signal  Estimates 

Disturbance  Potential  Estimate: 


[f]  =  [SEA]([cw]  0  (^mnlDi-SEFJ3, 

PxQ  PxK  KxL  KxL  LxQ 

Gravity  Vector  Estimates: 

[f]  =  [CLx]([s„n]  ®  [«„.,j)(ssy]T 

PxQ  PxK  KxL  KxL  LxQ 

[|]  =  |ssx]([am„]  ®  MJiesy]1' 

PxQ  PxK  KxL  KxL  LxQ 

f  ]  =  [5SA]([amn]  0  [0mn4])[ssy]r 

PxQ  PxK  KxL  KxL  LxQ 

Gravity  Gradient  Estimates: 

g]  =  [SSX]([Smn]  0  [0mn5])[SElT 

PxQ  PxK  KxL  KxL  LxQ 

g]  =  [S2A1([Sm„]®[«m„6])[S2yF 

PxQ  PxK  KxL  KxL  LxQ 

[§]  =  (s,sA'i([a„„i  ®  [«m„7])issyiT 

PxQ  PxK  KxL  KxL  LxQ 

[S]  =  (esA'l([s»»]  ®  MJlcsyf 

PxQ  PxK  KxL  KxL  LxQ 


(71) 

(72) 

(73) 

(74) 

(75) 

(76) 

(77) 

(78) 
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[Ufc]  =  (ss-WU  ®  [«„„9])|CEyf 

PxQ  PxK  KxL  KxL  LxQ 

[H]  =-  (CEX]([S„.)  ®  |«„„,ol)tSElT 

PxQ  PxK  KxL  KxL  LxQ 

where  <g>  has  been  defined  earlier  as  point-by-point  matrix  multiplication  and 


YX 

[5SX(i, ;)]  =  sin7r(t  x  j)-j- 

l  <»  <P 

1  <j<I< 

(81) 

ry 

|SSK(i,j)]T=8inir(ixj)^- 

l<i<L 

i  <j<Q 

(82) 

YX 

[CZX(i,j)]  =  cos  7r(t  x  j)— 

l<i<P 

1  <j<  K 

(S3) 

YY 

[CSK(i,i)]T  =  cosir(t  x  j)=£- 

1  <i<L 

1  <J  <Q 

(84) 

ft  .  —  ^  ~-Cmn  IM-01 

mnl  Vab 

(85) 

Bn  —  ^  n  g“"c»nn|^+-D| 

mn2  SAB 

(86) 

Q  ,  _  2  U  „-Cmn|/l+D| 

~  Sab  " 

(87) 

„  _  -2sgn (h  +  D)_  __ 

”mn4  •  \JAlB 

Cmn  |h-f  D| 

(ss) 

Be  —  ~  a? 

mnS  Sab  m 

(89) 

ft  „  —  — 2  A2„-emn|M-£>| 

mne  SABn 

(90) 

ft  _  -  2  c2  e-Cmn|/l+D| 

(91) 

0mn8  =  sfjQCLmhne  Cm"|,l+D| 

(92) 

_  -2sgn(fc  +  D)t  _ 

"mn9  \/  A.  B  °nCmn 

g—Cmn  |h4*£| 

(93) 

_  — 2sgn(/i +  D) 

“mnlO  —  \f~AB 

g— Ctnri 

(94) 

(95) 
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B  Software  Implementation 

B.l  Estimator  Module 

The  two  software  modules  KLE  (Karhunen-Loeve  estimator)  and  SYN  (synthesizer)  form 
the  two  stages  of  the  estimation  process.  The  estimator  equations  of  the  previous  section 
have  been  modified  in  form  only,  not  in  substance,  and  adapted  for  programming. 

In  the  first  stage  of  the  estimator,  KLE,  the  *th  measurement  field  is  sampled  on  the 
rectangular  grid.  It  is  transformed  on  the  left  and  right  by  the  sine  and/or  cosine  matrix 


given  by 

7 rkl 

M=i,-. 

.,K 

(96) 

2smI<  +  V 

or 

irkl 

M  =  i,. 

,.,I< 

(97) 

COoI<  +  V 

or 

irkl 

2sini+ 1- 

M  =  i,.. 

,.,L 

(98) 

or 

irkl 

2cosi+ 1’ 

..,.L 

(99) 

depending  upon  the  matrix  type  as  given  in  table  7  and  the  number  of  grid  samples  in  the 


Table  7:  Signal  dependent  functions 


ith  signal 

Ikli 

left 

right 

T 

2 

sine 

sine 

Tx 

2ak 

cosine 

sine 

Ty 

26/ 

sine 

cosine 

Tz 

-2  cm 

sine 

sine 

TXx 

-M 

sine 

sine 

Tyy 

-26? 

sine 

sine 

Tzz 

2  cl, 

sine 

sine 

Txy 

2akbi 

cosine 

cosine 

T» 

-26/q/ 

sine 

Txz 

—2  (ikCki- 

sine 

interior  of  the  survey  region,  K  in  the  x  direction  and  L  in  the  y  direction.  The  multipler  of 
2  is  for  compatibility  with  the  FST  and  FCT  functions  in  the  Swarztrauber  FFT  software 
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package  FFTPACK.  If  the  table  entry  is  “cosine”,  then  the  cosine  transform  is  followed  by 
an  application  of  the  matrix 


6ij  +  2|  i  — 

j  +  lb,  ij  =  1 

(100) 

6{j  +  2| i  - 

*»  j  —  1)  •••,  L 

(101) 

Note  that  |fc|2  denotes  fc  modulo  2,  i.e.  1  for  k  odd  and  0  for  k  even  and  8  is  the  Kronecker 
delta.  The  transformed  ith  measurement  field  is  denoted  by  Zku. 


The  following  equations  combine  the  transformed  measurement  fields  to  yield  the  K-L 
coefficients  a*/.  _ 


Cki  =  \]a\  +  ak  =  7 rk/A,  bt  =  irl/B 

(102) 

tki  =  exp  (-cw  •  ( h  +  D))\/~AB 

(103) 

hi  =  cro/cki° 

(104) 

o  _  hi  •  (tkijkii) 

[l  +  4f  Aw  Y.j(tkiykij)2/<r]\ 

(105) 

AB/ 4  ^ 

akl-  (K  +  l)(L  +  l)X'/3khZkh 

(106) 

The  main  computational  loop  in  the  program  increments  the  index 
merits  are  processed  sequentially.  Only  the  selected  measurements  are 
coefficient  field  aki  is  stored  in  a  disk  file. 

i  so  that  measure- 
processed.  The  K-L 

B.2  Synthesizer  Module 

The  SYN  module  reads  the  K-L  coefficient  field  otkh  recomputes  tki  and  7 ku  and  computes 
the  field 

z'ku  =  tkVykii<Xkt/4 

(107) 

where  now  i  represents  the  quantity  to  be  estimated.  The  z'kli  matrix  for  the  given  i  is  then 
multiplied  on  the  left  and  right  by  the  appropriate  sine  and/or  cosine  matrix  to  produce 
the  estimate  of  the  ith  signal. 


B.3  Matrix  Formulation 

It  is  useful  to  look  at  the  computational  process  in  high  level  matrix  notation.  In  fact,  the 
computations  can  be  summarized  in  just  two  equations: 

A  =  'EBiO{ViZiR'i)  (IQS) 
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and 

Yi  =  L,  {A  ©  r,-}  R,  (109) 

In  equation  (108)  the  *th  measurement  field  is  transformed  on  the  left  and  the  right  by 
L'f  and  R',-  which  are  the  sine  or  cosine  matrix  of  Equations  (96)-(99)  multiplied  by  the 
matrix  of  Equations  (100)-(10l).  This  puts  the  measurement  field  into  the  K-L  domain 
where  it  is  multiplied  by  the  weighting  matrix  The  elements  of  I?,-,  which  may  be 
inferred  from  Equation  (106)  and  its  predecessors,  are  given  by 


IA)h  = 


AB/4  ,, 
(A’  +  l)(£  +  l)' 


(110) 


The  “©”  product  is  simply  the  element  by  element  product  given  by 


A  ©  B  =  C  akibki  =  cm 


(111) 


The  summation  over  the  measurement  index  i  yields  the  field  of  K-L  coefficients.  In 
Equation  (109)  we  obtain  an  estimate  of  the  ith  signal  field.  The  K-L  coefficient  matrix 
is  multiplied  element  by  element  by  another  weighting  matrix  I\,  which  may  be  inferred 
from  Equation  (107)  and  its  predecessors.  Specifically,  the  elements  of  I\  are  given  by 


(112) 


The  result  is  left  and  right  multiplied  by  L;  and  R,-,  respectively,  each  of  which  is  either 
the  sine  or  cosine  matrix  of  Equations  (96)-(99). 

It  is  the  power  of  the  Karhunen-Loeve  method  that  the  estimation  equations  can  be 
cast  in  such  a  simple  form  as  Equations  (108)  and  (109),  in  particular  the  fact  that  the 
elements  of  the  K-L  coefficient  matrix  are  uncorrelated. 
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Figure  16:  Measured  Txy,  Single  DipoL 
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Figure  17:  Measured  Tyz,  Single  Dipole 
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Figure  18:  Measured  Txz ,  Single  Dipole 
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Figure20:  Tx  given  Txx,  Single  Dipole 
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Figure21:  Tv  given  Tjr,  Single  Dipole 
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Figure23:  a  given  TVJ/,  Single  Dipole 
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Figure28:  Tx  given  Tzz,  Single  Dipole 
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Figure30:  Tz  given  Tzz,  Single  Dipole 
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Figure32:  Tx  given  Txv,  Single  Dipole 
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Figure33:  T„  given  Txy,  Single  Dipole 
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Figure35:  &  given  Tyz\  Single  Dipole 
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Figure36:  Tx  given  I’m.  Single  Dipole 
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Figure38:  Tz  given  Tyi,  Single  Dipole 
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Figure39:  a  given  Txz,  Single  Dipole 
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Figure40:  T*  given  Txz,  Single  Dipole 
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Figure41:  fy  given  T„z,  Single  Dipole 
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Figure42:  Tz  given  Txz,  Single  Dipole 
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Figure43:  a  given  full  tensor  gradient,  Single  Dipole 


Tx  hat  (m/sec/sec) 


o 

II 

jC 

in 

<b 


E 

OJ 

II 

jC 

♦-> 

03 

V, 

O 

in 

c 

C b 


\ s 
o; 
a 
s 


a. 

c 

c 

r~ 

cz 

(D 

r— q 

CT 

c 

r— < 

co 


70 


Figure44:  Tx  given  full  tensor  gradient,  Single  Dipole 
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Figure45:  Tv  given  full  tensor  gradient,  Single  Dipole 
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Figure46:  Tz  given  full  tensor  gradient,  Single  Dipole* 
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Figure54:  Tx  given  T«,  File  1 
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Figure55:  Ty  given  Tir,  File  1 
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Figure56:  Tz  given  T„,  File  1 
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Figure57:  a  given  Tvy,  File  1 
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Figure58:  Tx  given  Tvy,  File  1 
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Figure59:  Ty  given  Tvv,  File  1 
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Figure60:  Tz  given  Tvv,  File  1 
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Figure61:  a  given  T«,  File  1 
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FigureC2:  Tx  given  T«,  File  1 
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Figure63:  Ty  given  T„,  File  1 
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FigureG4:  Tz  given  T«,  File  1 
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Figure65:  a  given  Txu,  File  1 
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Figure66:  Tx  given  TXV1  File  1 
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Figure67:  Ty  given  Txy,  File  1 
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Figure69:  a  given  TyI,  File  1 
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Figure70:  Tx  given  Tyz,  File  1 
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Figure71:  Tv  given  Tyz,  File  1 
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Figure72:  Tz  given  Tys,  File  1 
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Figure73:  a  given  Txz,  File  1 
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Figure74:  Tx  given  Txz,  File  1 
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Figure76:  Tz  given  Txz,  File  1 
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Figure7S:  Tx  given  full  tensor  gradient.  File  1 


Tg  hat  Cmgal) 


i - 1 — — - 1 - 1 


T-t 

ru 

n 

o 

o 

T 

n 

in 

n 

05 

O' 

o 

& 

3 

o 

*— • 

o 

r- 

n 

in 

OJ 

o 

II 

£ 

*-> 

m 

cu 


e 

vD 

o 

II 

£ 

♦-> 

03 

\_ 

O 

W 

c 

QJ 

♦-> 


D 


in 

03 

Qj 

e 


Qj 


a j 

a 

03 

i- 


105 


Figure79:  Ty  grven  full  tensor  gradient,  File  1 
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Figure80:  Tz  given  full  tensor  gradient,  File  1 
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Figure84:  Measured  T-y,  File  2 
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Figure85:  Measured  Tyz,  File  2 
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FigureSG:  Measured  Tx~,  File  2 
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Figure87:  a  given  TXI,  File  2 
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Figure89:  Ty  given  Txx,  File  2 
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Figure90:  Tx  given  Txx,  File  2 
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Figure91:  a  given  Tyy,  File  2 
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Figure92:  Tr  given  Tyy,  File  2 
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Figure95:  a  given  T«,  File  2 
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FigureOG:  Tx  given  T«,  File  2 
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Figure97:  fy  given  T«,  File  2 
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FigureOS:  Tz  given  Tzz,  File  2 
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Figure99:  a  given  TX]n  File  2 
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FigurelOO:  fx  given  Txv,  File  2 
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FigurelOl:  Tv  given  Tiy,  File  2 


Figurel02:  Tz  given  Txv,  File  2 
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Figure  103:  a  given  TyZ)  File  2 
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Figurel05:  Ty  given  Tyt,  File  2 
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Figurel06:  Tz  given  Tvz,  File  2 
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Figurel07:  a  given  Txz,  File  2 
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Figurel08:  Tx  given  T«,  File  2 
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FigureliO:  Tz  given  Txz,  File  2 
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Figurelll:  a  given  full  tensor  gradient,  File  2 


Tx  hat  (mgal  ) 


in 

o 

in 

a 

o 

in 

o 

U) 

on 

in 

00 

*-4 

CO 

ir 

li- 

=f 

T 

CU 

ON 

r— * 

i-i 

CO 

I  I 


O 

II 

H 


<-> 

in 

CL) 


£ 

.X 

vD 

CD 

II 

r 

03 


0 

in 

c 

CD 


m 

(0 

<b 

£ 


cu 

Qj 


4- 


(D 

a 

0] 

i- 


138 


Figurell2:  Tx  given  full  tensor  gradient,  File  2 


Figurell3:  Ty  given  full  tensor  gradient,  File  2 
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Figurell4:  Tz  given  full  tensor  gradient,  File  2 
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Figurell5:  a  given  full  tensor  gradient.  Files  1  &  2,  intersection  method 


Tx  given  full  tensor  gradient,  Files  1  &  2,  intersection  method 
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Figurell7:  Tv  given  full  tensor  gradient.  Files  1  &  2,  intersection  method 
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Figurell8:  Tz  given  full  tensor  gradient,  Files  1  &  2,  intersection  method 
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Figurel20:  Tx  given  full  tensor  gradient,  Files  1  &  2,  union  method 


Ty  hat  Cmgal ) 


A 


147 


Tape .  1 .»  files  1  &  2  combined,  union  method 
Figurel21:  Ty  given  full  tensor  gradient.  Files  1  &  2,  union  method 
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Figurel22:  Tz  given  full  tensor  gradient,  Files  1  &  2,  union  method 
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Figurel23:  Truth  value  of  Tx,  File  3 
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Figurel24:  Truth  value  of  T„,  File  3 
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Figurel25:  Truth  value  of  Tz,  File  3 


